library(readxl) library(readODS) library(writexl) library(ggthemes) library(lvplot) library(glue) library(tidyverse) library(plotly) # gráfico 3d library(randtests) library(goft)
Tópico 5
library(readxl) library(readODS) library(writexl) library(ggthemes) library(lvplot) library(glue) library(tidyverse) library(plotly) # gráfico 3d library(randtests) library(goft)
O gerador congruente linear mais famoso é o randu proposto em 1951 e implementado pela IBM logo em seguida.
Restrições
Suponha que \(s_{k+1} < 2^{31}\), então \[ \begin{split} s_{k+2} &= (2^{16} + 3) s_{k+1} = (2^{16} + 3) (2^{16} + 3) s_{k},\\ &= (2^{32} + 6 \cdot 2^{16} + 9) s_k = 6 (2^{16} + 3) s_k - 9 s_k,\\ &= 6 s_{k+1} - 9 s_k. \end{split} \]
Ou seja, existe uma relaçãos simples entre \(x_{k+2}\), \(x_{k+1}\) e \(x_k\), e podemos descobrir o gerador pseudoaleatório a partir da sequência de número pseudoaleatórios.
randu <- function(n, seed = 1) {
x <- vector(mode = "double", n)
a <- 2^16 + 3
m <- 2^31
x[1] <- (a * seed) %% m
for (i in 2:n) x[i] <- (a * x[i - 1]) %% m
x / m
}
Marsaglia (1968) provou os valores desta matriz estão dentro de 15 hiperplanos.
Vamos usar o teste de aleatoriedade (Wald e Wolfowitz 1940): runs.test() do pacote randtests.
tamanho_amostra <- 9999 amostra <- randu(tamanho_amostra) matriz <- tibble( xk = amostra[seq_len(tamanho_amostra) %% 3 == 0], xk1 = amostra[seq_len(tamanho_amostra) %% 3 == 1], xk2 = amostra[seq_len(tamanho_amostra) %% 3 == 2] )
ks.test(amostra, "punif")
## ## One-sample Kolmogorov-Smirnov test ## ## data: amostra ## D = 0.0063368, p-value = 0.8168 ## alternative hypothesis: two-sided
runs.test(amostra)
## ## Runs Test ## ## data: amostra ## statistic = -2.0203, runs = 4899, n1 = 4999, n2 = 4999, n = 9998, p-value = 0.04335 ## alternative hypothesis: nonrandomness
plot_ly(matriz, x = ~xk, y = ~xk1, z = ~xk2, size = 0.2)
Definição 4 (Gerador Defasado de Fibonacci) Seja \(i,j,k\) números inteiros tal que \(0 < k < j < i\), e considere o Gerador Congruente Linear \(\Xi\) com em que \(S = \mathbb{Z} \cap [0, 2^{31}]\), \(T(s_i) = (s_{i-j} + s_{j-k}) \mod 2^{31}\), \(U = \frac{S}{2^{31}}\), \(G(s) = \frac{s}{2^{31}}\) e \(s_0\) é um conjunto de valores iniciais tal que \(\lvert s_0 \rvert > j\). Chamamos \(\Xi\) de Gerador Desafado de Fibonacci.
Este algoritmo é uma melhoria sobre o randu, e pode ser usado para simulações. (Mas não é recomendado para criptografia).
O R implementou o algoritmo Mersenne Twister (vide Härdle, Okhrin, e Okhrin 2017 para mais detalhes) com a função runif(n).
k <- 5
j <- 17
m <- 2^31
semente <- runif(j, 0, m) |> round()
gdf <- function(n, seed = semente) {
x <- seed
for (i in seq_len(n)) {
x[j + i] <- (x[i] + x[j + i - k]) %% m
}
x[(j + 1):length(x)] / m
}
tamanho_amostra <- 9999 amostra <- gdf(tamanho_amostra) matriz <- tibble( xk = amostra[seq_len(tamanho_amostra) %% 3 == 0], xk1 = amostra[seq_len(tamanho_amostra) %% 3 == 1], xk2 = amostra[seq_len(tamanho_amostra) %% 3 == 2] )
ks.test(amostra, "punif")
## ## One-sample Kolmogorov-Smirnov test ## ## data: amostra ## D = 0.0082316, p-value = 0.507 ## alternative hypothesis: two-sided
runs.test(amostra)
## ## Runs Test ## ## data: amostra ## statistic = 1.4002, runs = 5070, n1 = 4999, n2 = 4999, n = 9998, p-value = 0.1615 ## alternative hypothesis: nonrandomness
plot_ly(matriz, x = ~xk, y = ~xk1, z = ~xk2, size = 0.2)
ex_discreto <- function(n) {
if (length(n) > 1) stop("n precisa ser um valor escalar inteiro.")
seq_len(n) |>
map_dbl(\(x) {
u <- runif(1)
if(u <= 0.3) {
return(0)
} else if (0.3 < u && u <= 0.5) {
return(1)
} else if (u > 0.5) {
return(2)
}
})
}
set.seed(121343) tibble(amostra = ex_discreto(1000)) |> group_by(amostra) |> summarise(freq = n()) |> mutate(fr = freq / sum(freq))
## # A tibble: 3 × 3 ## amostra freq fr ## <dbl> <int> <dbl> ## 1 0 304 0.304 ## 2 1 195 0.195 ## 3 2 501 0.501
lambda <- 2
exp_2 <- function(n) {
if (length(n) > 1) stop("n precisa ser um valor escalar inteiro.")
seq_len(n) |>
map_dbl(\(x) {
u <- runif(1)
if (u <= 0) {
return(0)
} else {
return(-log(1 - u) / lambda)
}
})
}
n <- 1000 k <- (1 + log2(n)) |> round() dados <- tibble(amostra = exp_2(n)) ggplot(data = dados) + geom_histogram(mapping = aes(x = amostra, y = ..density..), bins = k) + theme_minimal() + labs(x = "Amostra", y = "Densidade de Frequência")
ks.test(dados$amostra, pexp, rate = 2)
## ## One-sample Kolmogorov-Smirnov test ## ## data: dados$amostra ## D = 0.03783, p-value = 0.1143 ## alternative hypothesis: two-sided
Geralmente, usamos \(c\) como sendo o valor máximo de \(\frac{f(x)}{g(x)}\).
beta_21 <- function(n) {
if (length(n) > 1) stop("'n' precisa ser um escalar inteiro.")
seq_len(n) |>
map_dbl(\(i) {
y <- runif(1); u <- runif(1)
while(u > dbeta(y, shape1 = 2, shape2 = 1) / (2 * dunif(y))) {
y <- runif(1); u <- runif(1)
}
return(y)
})
}
n <- 1000 k <- (1 + log2(n)) |> round() dados <- tibble(amostra = beta_21(n)) ggplot(data = dados) + geom_histogram(mapping = aes(x = amostra, y = ..density..), bins = k) + theme_minimal() + labs(x = "Beta(2, 1)", y = "Densidade de Frequência")
ks.test(dados$amostra, pbeta, shape1 = 2, shape2 = 1)
## ## One-sample Kolmogorov-Smirnov test ## ## data: dados$amostra ## D = 0.030908, p-value = 0.295 ## alternative hypothesis: two-sided
| \(x\) | 1 | 2 | 3 | 4 | 5 |
|---|---|---|---|---|---|
| \(f(x)\) | \(0,15\) | \(0,22\) | \(0,33\) | \(0,10\) | \(0,20\) |
fp <- c(0.15, 0.22, 0.33, 0.10, 0.20)
const <- 1.65
mar_discreto <- function(n) {
if (length(n) > 1) {
stop("'n' precisa um escalar inteiro.")
}
seq_len(n) |>
map_dbl(\(i) {
y <- sample.int(5, 1); u <- runif(1)
while(u > fp[y] / (0.2 * const)) {
y <- sample.int(5, 1); u <- runif(1)
}
return(y)
})
}
n <- 1000 dados <- tibble(amostra = mar_discreto(n)) print(fp)
## [1] 0.15 0.22 0.33 0.10 0.20
dados |> group_by(amostra) |> summarise(freq = n()) |> mutate(fr = freq / sum(freq))
## # A tibble: 5 × 3 ## amostra freq fr ## <dbl> <int> <dbl> ## 1 1 146 0.146 ## 2 2 211 0.211 ## 3 3 340 0.34 ## 4 4 98 0.098 ## 5 5 205 0.205
Sejam \(U_1\sim U(0,1)\) e \(U_2\sim U(0,1)\) duas variáveis aleatórias contínuas independentes, e considere \[ \begin{split} X_1 &= \sqrt{-2\ln(U_1)}\cos(2\pi U_2),\\ X_2 &= \sqrt{-2\ln(U_1)}\sin(2\pi U_2), \end{split} \] então \(X_1\) e \(X_2\) são independentes, e \(X_1 \sim N(0,1)\) e \(X_2 \sim N(0,1)\).
box_muller <- function(n) {
if (length(n) > 1) stop("'n' precisa ser um escalar inteiro.")
seq_len(n) |>
map_dbl(\(i) {
u1 <- runif(1); u2 <- runif(1)
sqrt(-2 * log(u1)) * cos(2 * pi * u2)
})
}
n <- 1000 k <- (1 + log2(n)) |> round() dados <- tibble(amostra = box_muller(n)) ggplot(data = dados) + geom_histogram(mapping = aes(x = amostra, y = ..density..), bins = k) + theme_minimal() + labs(x = "Distribuição normal.", y = "Densidade de Frequência")
ks.test(dados$amostra, pnorm)
## ## One-sample Kolmogorov-Smirnov test ## ## data: dados$amostra ## D = 0.022602, p-value = 0.6866 ## alternative hypothesis: two-sided
Sejam \(U_1\sim U(-1,1)\) e \(U_2\sim U(-1,1)\) duas variáveis aleatórias contínuas independentes. Suponha que \(\omega = U_1^2 + U_2^2 \leq 1\), e considere \[ \begin{split} X_1 &= \sqrt{\frac{-2\ln(\omega)}{\omega}} U_1,\\ X_2 &= \sqrt{\frac{-2\ln(\omega)}{\omega}} U_2. \end{split} \] Então é possível provar que \(X_1\) e \(X_2\) são duas variáveis aleatórias independentes, e \(X_1 \sim N(0, 1)\) e \(X_2 \sim N(0,1)\).
Para gerar valores de \(U(-1, 1)\) usamos o método da transformação inversa. Note que \(2 \cdot U -1 \sim U(-1, 1)\) com \(U \sim U(0,1)\).
metodo_polar <- function(n) {
if (length(n) > 1) stop("'n' precisa ser um escalar inteiro.")
seq_len(n) |>
map_dbl(\(i) {
u1 <- 2 * runif(1) - 1; u2 <- 2 * runif(1) - 1; w <- u1^2 + u2^2
while(w > 1) {
u1 <- 2 * runif(1) - 1; u2 <- 2 * runif(1) - 1; w <- u1^2 + u2^2
}
return(sqrt(-2 * log(w) / w) * u1)
})
}
n <- 1000 k <- (1 + log2(n)) |> round() dados <- tibble(amostra = metodo_polar(n)) ggplot(data = dados) + geom_histogram(mapping = aes(x = amostra, y = ..density..), bins = k) + theme_minimal() + labs(x = "Valores amostrados", y = "Densidade de Frequência")
ks.test(dados$amostra, pnorm)
## ## One-sample Kolmogorov-Smirnov test ## ## data: dados$amostra ## D = 0.021913, p-value = 0.7229 ## alternative hypothesis: two-sided
gamma_alpha_grande <- function(n, shape, scale) {
if (length(n) > 1) stop("'n' precisa ser um escalar inteiro.")
if (shape <= 1) stop("'shape' precisa ser um escalar maior que 1.")
if (scale < 0) stop("'scale' precisa ser um escalar positivo.")
seq_len(n) |>
map_dbl(\(i) {
a <- 1 / sqrt(2 * shape - 1); b <- shape - log(4)
ce <- shape + 1 / a
u1 <- runif(1); u2 <- runif(1)
v <- a * log(u1 / (1 - u1)); x <- shape * exp(v)
while(b + ce * v - x < log(u1^2 * u2)) {
u1 <- runif(1); u2 <- runif(1)
v <- a * log(u1 / (1 - u1)); x <- shape * exp(v)
}
return(scale * x)
})
}
n <- 1000 k <- (1 + log2(n)) |> round() dados <- tibble(amostra = gamma_alpha_grande(n, shape = 2, scale = 4)) ggplot(data = dados) + geom_histogram(mapping = aes(x = amostra, y = ..density..), bins = k) + theme_minimal() + labs(x = "valores amostrados", y = "Densidade de Frequência")
ks.test(dados$amostra, pgamma, shape = 2, scale = 4)
## ## One-sample Kolmogorov-Smirnov test ## ## data: dados$amostra ## D = 0.020675, p-value = 0.7861 ## alternative hypothesis: two-sided
gamma_alpha_pequeno <- function(n, shape, scale) {
seq_len(n) |>
map_dbl(\(i) {
b <- (shape + exp(1)) / exp(1); u1 <- runif(1); u2 <- runif(1); p <- b * u1; E <- -log(u2)
laco <- TRUE
while(laco) {
if (p <= 1) {
x <- p^(1 / shape)
if (x > E) {
b <- (shape + exp(1)) / exp(1); u1 <- runif(1); u2 <- runif(1); p <- b * u1; E <- -log(u2)
} else {
laco <- FALSE
}
} else {
x <- -log((b - p) / shape)
if ((1 - shape) * log(x) > E) {
b <- (shape + exp(1)) / exp(1); u1 <- runif(1); u2 <- runif(1); p <- b * u1; E <- -log(u2)
} else {
laco <- FALSE
}
}
}
return(scale * x)
})}
n <- 1000 k <- (1 + log2(n)) |> round() dados <- tibble(amostra = gamma_alpha_pequeno(1000, shape = 0.5, scale = 4)) ggplot(data = dados) + geom_histogram(mapping = aes(x = amostra, y = ..density..), bins = k) + theme_minimal() + labs(x = "valores amostrados", y = "Densidade de Frequência")
ks.test(dados$amostra, pgamma, shape = 0.5, scale = 4)
## ## One-sample Kolmogorov-Smirnov test ## ## data: dados$amostra ## D = 0.024642, p-value = 0.5782 ## alternative hypothesis: two-sided
Seja \(X \sim B(\alpha, \beta)\) com função densidade dada por \(f(x) = \frac{\Gamma(\alpha) \Gamma(b)}{\Gamma(\alpha + \beta)} x^{\alpha - 1} (1 - x)^{\beta -1}\).
r_beta <- function(n, shape1, shape2) {
seq_len(n) |>
map_dbl(\(i) {
u1 <- runif(1); u2 <- runif(1); v1 <- u1^(1 / shape1)
v2 <- u2^(1 / shape2); w <- v1 + v2
while(w > 1) {
u1 <- runif(1); u2 <- runif(1); v1 <- u1^(1 / shape1)
v2 <- u2^(1 / shape2); w <- v1 + v2
}
return(v1 / w)
})
}
n <- 1000 k <- (1 + log2(n)) |> round() dados <- tibble(amostra = r_beta(n, shape1 = 2, shape2 = 3)) ggplot(data = dados) + geom_histogram(mapping = aes(x = amostra, y = ..density..), bins = k) + theme_minimal() + labs(x = "valores amostrados", y = "Densidade de Frequência")
ks.test(dados$amostra, pbeta, shape1 = 2, shape2 = 3)
## ## One-sample Kolmogorov-Smirnov test ## ## data: dados$amostra ## D = 0.035974, p-value = 0.1502 ## alternative hypothesis: two-sided
Seja \(X \sim b(n, p)\).
r_binom <- function(n, size, p) {
seq_len(n) |>
map_dbl(\(i) {
seq_len(size) |> map_dbl(~ runif(1) <= p) |> sum()
})
}
n <- 1000 k <- (1 + log2(n)) |> round() dados <- tibble(amostra = r_binom(1000, 5, 0.75)) ggplot(data = dados) + geom_bar(aes(x = amostra, y = 100 * ..prop..)) + theme_minimal() + labs(x = "valores amostrados", y = "Porcentagem")
n <- 1000 size <- 5; p <- 0.55 suporte <- 0:size fda <- stepfun(x = suporte, y = c(0, pbinom(suporte, size = size, prob = p))) dgof::ks.test(rbinom(1000, size, p), fda)
## ## One-sample Kolmogorov-Smirnov test ## ## data: rbinom(1000, size, p) ## D = 0.024218, p-value = 0.6006 ## alternative hypothesis: two-sided
Seja \(X \sim Poison(\lambda)\).
r_pois <- function(n, lambda) {
seq_len(n) |>
map_dbl(\(i) {
k <- 0; p <- 1
while(exp(-lambda) < p) {
k <- k + 1; p <- p * runif(1)
}
return(k - 1)
})
}
n <- 1000 k <- (1 + log2(n)) |> round() dados <- tibble(amostra = r_pois(1000, lambda = 4)) ggplot(data = dados) + geom_bar(aes(x = amostra, y = 100 * ..prop..)) + theme_minimal() + labs(x = "valores amostrados", y = "Porcentagem")
suporte <- 0:10000 fda <- stepfun(suporte, c(0, ppois(suporte, lambda = 4))) dgof::ks.test(dados$amostra, fda)
## ## One-sample Kolmogorov-Smirnov test ## ## data: dados$amostra ## D = 0.026163, p-value = 0.5004 ## alternative hypothesis: two-sided
r_mnorm <- function(n, sigma, media) {
media <- matrix(media, ncol = 1)
R <- chol(m_var)
seq_len(n) |>
sapply(\(i){
return(R %*% matrix(nrow(sigma) |> rnorm(), ncol = 1) + media)
}) |>
t()
}
p <- 0.90
dim <- 10
m_var <- seq_len(dim) |>
sapply(\(i) {
seq_len(dim) |> map_dbl(\(j) p^abs(i - j))
})
v_media <- rep(2, dim)
dados <- r_mnorm(1000, sigma = m_var, media = v_media)
goft::mvshapiro_test(dados)
## ## Generalized Shapiro-Wilk test for Multivariate Normality ## ## data: dados ## MVW = 0.99832, p-value = 0.5471
Ahrens, Joachim H, e Ulrich Dieter. 1974. «Computer methods for sampling from gamma, beta, poisson and bionomial distributions». Computing 12 (3): 223–46.
Atkinson, AC, e MC Pearce. 1976. «The computer generation of beta, gamma and normal random variables». Journal of the Royal Statistical Society: Series A (General) 139 (4): 431–48.
Cheng, RCH. 1977. «The generation of gamma variables with non-integral shape parameter». Journal of the Royal Statistical Society: Series C (Applied Statistics) 26 (1): 71–75.
Härdle, Wolfgang Karl, Ostap Okhrin, e Yarema Okhrin. 2017. Basic elements of computational statistics. Springer.
Marsaglia, George. 1968. «Random numbers fall mainly in the planes». Proceedings of the National Academy of Sciences of the United States of America 61 (1): 25.
Neave, Henry R. 1973. «On using the Box-M üller transformation with multiplicative congruential pseudo-random number generators». Journal of the Royal Statistical Society: Series C (Applied Statistics) 22 (1): 92–97.
Wald, Abraham, e Jacob Wolfowitz. 1940. «On a test whether two samples are from the same population». The Annals of Mathematical Statistics 11 (2): 147–62.